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> Abstract 

We present the results of high-statistics calculations of correlation functions generated with single- 
baryon interpolating operators on an ensemble of dynamical anisotropic gauge-field configurations 
generated by the Hadron Spectrum Collaboration using a tadpole-improved clover fermion action 
and Symanzik-improved gauge action. A total of 292,500 sets of measurements are made using 
1194 gauge configurations of size 20^ x 128 with an anisotropy parameter = hs/ht = 3.5, a spatial 
lattice spacing of hg = 0.1227 it 0.0008 fm, and pion mass of M^^ ~ 390 MeV. Ground state baryon 
>■ masses are extracted with fully quantified uncertainties that are at or below the ~ 0.2%-level in 

^ lattice units. The lowest-lying negative-parity states are also extracted albeit with a somewhat 

H lower level of precision. In the case of the nucleon, this negative-parity state is above the Ntt 

threshold and, therefore, the isospin-2 vrA^ s-wave scattering phase-shift can be extracted using 
Liischer's method. The disconnected contributions to this process are included indirectly in the 
gauge-field configurations and do not require additional calculations. The signal-to-noise ratio in 
the various correlation functions is explored and is found to degrade exponentially faster than 
naive expectations on many time-slices. This is due to backward propagating states arising from 
the anti-periodic boundary conditions imposed on the quark-propagators in the time-direction. 
We explore how best to distribute computational resources between configuration generation and 
propagator measurements in order to optimize the extraction of single baryon observables. 
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I. INTRODUCTION 



One of the primary goals of lattice QCD (LQCD) is to calculate the properties and inter- 
actions of nucleons and, more generally, systems comprised of multiple hadrons. Precise 
exploration of the simplest multi-hadron systems has recently become possible with sig- 
nificant advances in computing resources, as well as through algorithmic and theoretical 
developments. The two-pion system 7r+7r+ is the simplest of such multi-hadron systems to 
calculate in LQCD, and current computational resources have allowed for a precise deter- 
mination of the 7r+7r+ scattering length [Tl|2] at the ~ 1% level. Recently, we have explored 
systems comprised of up to twelve vr+'s [21 H] and also systems comprised of up to twelve 
K~^^s [5] for the first time, allowing a determination of the three- tt"*" and three- if ^ interac- 
tions. In general, a determination of the two-particle scattering amplitude, or multi-body 
interactions, with LQCD requires calculating the energy-eigenvalues of the system in the 
finite- volume [SI El El E] • The energy differences between the multi-particle energy-levels in 
the finite-volume and the sum of the particle masses determines the scattering amplitude or 
interaction. As processes of interest to low-energy nuclear physics are in the MeV energy- 
regime, while the masses of the baryons and nuclei are in the GeV regime, the energy-levels 
in the volume must be determined to high precision to yield useful constraints and pre- 
dictions for scattering amphtudes, phase-shifts and electroweak properties. Consequently, 
correlation functions of systems comprised of more than one hadron must be calculated with 
small statistical and systematic uncertainties (<^ 1 %) in order to provide useful information 
about low-energy nuclear interactions and nuclei. 

The correlation functions associated with systems of baryons (and, more generally, states 
other than the pion) suffer from an exponential degradation of the signal-to-noise ratio as 
a function of time as argued by Lepage [10]. The scale that dictates this degradation is 
the difference between the total energy of the baryons in the system and half of the total 
energy of hadrons that contribute to the correlation function associated with the square of 
the interpolating operator for the baryon system. An example is provided by the two point 
nucleon correlator (A^ is an interpolating field with the quantum numbers of the nucleon) 
where, 

signal ~ {N{t)N{0)) *^ Ze"^^*, (1) 
noise ~ ^J{N{t) N{t)N{0) iV(0)) *=5 Z'e-t^^-*, 

neglecting effects of the finite temporal extent which we discuss below (here Z and Z' are 
overlap factors). Since — fM^r > 0, the signal degrades exponentially in time with this 
exponent. Further, for multi-baryon systems, this exponent is scaled by the baryon number, 
B, and it consequently requires exponentially larger computational resources to calculate 
the properties of systems containing B > 1 baryons than a single baryon. In many regards, 
it is this signal-to-noise problem that distinguishes LQCD calculations of quantities typically 
of importance to nuclear physics from those typically of importance to particle physics. 

The main motivation for our present work is to explore very high statistics calculations of 
the energy spectrum of -B = 0, 1 correlation functions, quantifying the statistical scaling and 
identifying any issues that appear in the regime of precision calculations. More generally, we 
aim to assess the feasibility of extracting precise phase-shifts and multi-nucleon interactions 
from multi-baryon systems but we leave these discussions to subsequent work. Our focus is 
on the statistical scaling behavior of these measurements instead of on measuring physically 
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relevant quantities. Consequently, we work with a single ensemble of gauge configurations 
that was produced by the Hadron Spectrum Collaboration [TT] (the details of these config- 
urations are discussed below). The analysis presented here enables us to identify a number 
of issues that will be important to LQCD calculations of quantities where exponentially 
degrading signal-to-noise ratios dominant concern: 

1. While the classic argument of Lepage [T0| concerning the behavior of the signal-to-noise 
ratios of baryon correlation functions seems to be on a solid theoretical footing, it has 
yet to be explored and verified through direct calculation. We examine the signal-to- 
noise ratios of the single hadron correlation functions in detail and present a modified 
version of the Lepage argument that incorporates the finite extent of the temporal 
direction of the gauge-field configuration, focusing on the case of quark propagators 
subject to anti-periodic temporal boundary conditions (BCs). Over large regions of the 
temporal extent of the lattice, the signal-to-noise ratio degrades exponentially faster 



than expected from the original Lepage argument, see Sec. VIII 



2. At present, and even more so in the past, the generation of gauge-field ensembles 
consumes most of the computational resources of LQCD calculations.^ However, it 
is not clear what the optimal distribution of computational resources between gauge- 
field production and measurements (propagator calculations and contractions) is when 
one is interested in noisy quantities. To address this, we explore what can be accom- 
plished by performing hundreds of measurements per configuration, and how precisely 
the baryon ground-state masses can be determined from an ensemble of 1194 configura- 
tions. We also study whether the measurements "saturate" after some critical number 
have been performed on one configuration (that is, exhibit little or no improvement in 
uncertainties after a critical number of measurements), finding for baryons that they 
do not, even up to ~ 200 measurements per configuration. 

3. Correlation functions that are determined to high precision are amenable to analysis 
with a variety of techniques, beyond those typically used successfully with low statistics 
data. On these anisotropic configurations, multiple (five or more) exponential fits 
to such correlation functions become stable as statistical fluctuations decrease, and 
the ground state energies can be extracted with high precision. We show that the 
generalized effective mass (EM) method, in which multiple energies are extracted 
from a linear system (a method developed by Gaspard Riche de Prony in 1795) also 
becomes useful for correlation functions with small uncertainties. As two (different 
but correlated) correlation functions are computed per species of hadron, this method 
is extended to construct the matrix-Prony method, which is found to be a very clean 
and effective tool for determining the ground-state energies. 

4. While the correlation function generated by a single baryon interpolating operator will 
be dominated by the baryon ground-state at large times, it also contains contributions 
from all states that can couple to the operator. This includes multi-hadron states. The 
backward propagating component of the nucleon correlation function is dominated by 
the lowest energy negative-parity I = ^ state for the projection-operator we have 
applied to the correlation functions. By measuring the energy of this state, which is 



^ For example, the USQCD collaboration used ~ 60% of its resources for ensemble generation in 2008/9. 
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above the vrA^ threshold and therefore is a continuum state, the phase-shift associated 
with the s-wave vrA^ interaction is determined at this energy. The important point here 
is that this process contains disconnected diagrams, which are encoded in the gauge- 
field configurations, and do not require additional (of order the volume in number) 
calculations. 

5. We also note that thermal states, while strongly suppressed, are seen in our high pre- 
cision data. In these states, some part of the hadronic state propagates backward in 
time and can consequently manifest itself in the correlation function as an exponential 
with energy less than that of the zero temperature ground state. These contributions 
have amplitudes that are exponentially suppressed by the temporal extent of the con- 
figuration, but they can be extracted in certain temporal regions of the correlation 
function(s) where other components are also small. They can lead to pollution of the 
ground state signal. 

The structure of this work is as follows. Section HIl introduces the details of the lattice 



calculations we perform, and in section III we discuss our expectations for the hadron spec- 



trum on this ensemble. Section IV introduces the tools used in our analysis and presents 



detailed comparisons of the different methods we utilize. Following this, sections VI 



and VII present our main results for the pseudo-scalar mesons, ground-state baryons and 



negative-parity excited states, respectively. In sections VIII and IX we discuss the behaviour 
of noise in our measurements and investigate the scaling of uncertainties in hadron masses 
for varying numbers of gauge configurations and measurements. We conclude in section |X} 
In subsequent works, we will address states with baryon number, B > 1. 



II. LATTICE QCD CALCULATIONS 



A. Calculational Details 



In this study, we employ an ensemble of the = 2 + 1-fiavor anisotropic clover gauge-field 
configurations that are currently being produced by the Hadron Spectrum Collaboration [TT] . 
These ensembles are being generated with dynamical tadpole-improved clover fermions and 
a Symanzik-improved gauge action (see Ref [IT] for full details). All of the calculations 
that we present here were performed on a single ensemble of gauge-field configurations of 
size 20^ X 128 with an anisotropy parameter of ^ = bg/bt = 3.5, a spatial lattice spacing of 
h = bs = 0.1227 ± 0.0008 fm, a pion mass of ~ 390 MeV and a kaon mass of Mk ~ 
546 MeV. The ensemble used in this study contains 1194 configurations taken at intervals 
of 10 trajectories, after allowing 1000 trajectories for thermalization. Ref. [TT] provides a 
comprehensive analysis of autocorrelation times and thermalization. Some correlation is 
seen to be present between configurations separated by 30 trajectories. 

The light and strange quark-propagators were computed using the same fermion action 
used in the gauge-field generation. We use the clover discretization of the fermion action 
as it requires significantly less computational resources than, for instance, the Domain- Wall 
discretization, in both the production of gauge-field configurations and in the calculation 
of quark-propagators, while retaining an (9(6)-improved spectrum. Unlike the Domain- 
Wall discretization, the Clover discretization does not have a lattice chiral symmetry. At 
moderate lattice spacings, this may significantly impact the extraction of the properties and 
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FIG. 1: The number of propagators, A'prop, used in measurements of correlation functions on 
each configuration used in this study. For the purpose of clarity, bins of four configurations (40 
trajectories) have been averaged. The ensemble-average of the number of propagators calculated 
per configuration, 245, is indicated by the dashed horizontal line. 

interactions of pions and kaons, but it is not expected to produce systematic uncertainties 
that are as significant in the properties and interactions of baryons(this remains to be verified 
and will not be addressed here). 

The quark propagators are calculated with anti-periodic BC's imposed on the time- 
direction and periodic BCs imposed on the spatial directions. As multiple propagators 
are calculated on each configuration, iterative solvers beyond the simple conjugate gradient 
algorithm can provide significant speed improvements. In particular, we employ the deflated 
conjugate gradient algorithm proposed in Ref. [12], and implemented in the Chroma lattice 
field theory library [13j as the EigCG inverter. In our typical production runs, we compute 
from 30 to 100 propagators in sequence, observing a factor of ~ 7 improvement in inversion 
speed after the first few solves are used to deflate low-lying eigenvalues from subsequent 
inversions. Figure [T] shows the details of the propagators computed in this work. The his- 
togram indicates the number of propagators computed on each of the 1194 configurations 
(averaged over four adjacent configurations for clarity) . The total number of propagators 
computed in this dataset is 292,500, an average of 245 propagators per configuration (we 
note that the maximum number of point-to-all propagators that could be computed on each 
of the configurations is ~ 10^) 

Each propagator is generated from a gauge-invariantly Gaussian-smeared source [m [15] , 
on a stout-smeared gauge-field in order to optimize the overlap onto the ground-state 
hadrons. On each configuration, the locations of the propagator source points are chosen 
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FIG. 2: The separation between pairs of sources on a given configuration, defined to be the mini- 
mum distance between two sources, including the effect of the (anti-) periodic boundary conditions. 
The height of the bar at ii = corresponds to the total number of propagators. 



randomly throughout the configuration. In fig. |2| we show a histogram of the 4d-separation, 
R, between the each pair of sources on each configuration. The shoulder at i? ~ 4 fm ap- 
pears because of the (anti) periodic boundary conditions. The average source separation is 
(R) ~ 2.9 fm and the source density is 3.43 fm~^. 



B. Correlation Functions 

The propagators are used to compute two-point correlation functions which, for baryons, 
take the form 

ChAp; t) = Yl ^'"^^ < ^)^-(xo, 0) ) , (2) 

X 

where 7i°(x, t) is an interpolating operator for the appropriate baryon state, e.g., for the 
proton 7i"(x, t) = eabc {u°"^ C 75^'') m'^'"(x, t) where C is the charge conjugation matrix. The 
Dirac matrix F is an arbitrary particle-spin-projector and the point xq is the propagator 
source point. Similar correlation functions are used for the mesons. The interpolating- 
operator at the source, TC, is constructed from gauge-invariantly-smeared quark field opera- 
tors, while at the sink, the interpolating operator is constructed from either local quark field 
operators, or from the same smeared quark field operators used at the source, leading to two 
sets of correlation functions. For brevity, we refer to the two sets of correlation functions 
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that result from these source and sink operators as smeared-point (SP) and smeared- smeared 
(SS) correlation functions, respectively. 

We calculate the smeared-point and smeared-smeared correlation functions associated 
with the 7r+, ( J'^ = 0~) mesons, and the N, A, E, S ( J'^ = baryons. For the baryons, 
the energy projectors T± = ^{1 ±70) are used to project separately onto either the positive- 
or negative-energy (parity) states (one of which can be time-reversed and added to the other 
to improve statistics). Correlation functions associated with a given pair of interpolating 
fields are averaged over all sources on each configuration, producing one correlation function 
per interpolating operator pair per configuration. 

C. Statistical Behavior 

Before extracting results for observables, we analyze the statistical behavior of the mea- 
sured correlators. As the computational cost of each measurement is much less than the 
computational cost of generating each configuration, performing multiple, (9(10), measure- 
ments on each configuration is a practical way to cheaply reduce uncertainties and is an 
approach that has been used by many groups. Averaging the measurements on a given 
configuration produces a more accurate estimation of the correlation function on that con- 
figuration. A priori, one might argue that performing a significantly larger number, say 
(9(100-1000), of measurements on a given configuration is an inefficient use of computing 
resources as the additional measurements will contain little or no new information and will 
not decrease the statistical uncertainty in the measurements of interest. This argument is 
likely true for configurations extending over small volumes. Physically, one expects there 
are a number of length scales associated with the possible "saturation" density of the mea- 
surements on a given configuration. As the lightest hadron is the pion, one expects the 
critical saturation density of measurements to depend parametrically upon the dimension- 
less quantity p/M^, where p = Ngrc/V where V is the four-volume, and Nsrc is the number 
of measurements on the configuration. For a simple quantity such as the energy, E, of an 
eigenstate in the volume, one also expects to find a dependence upon p/E'^. For instance, 
we expect a dependence upon the dimensionless quantity A^src/(^ M^) for the nucleon. 
Figure [3] shows the scaling of the uncertainty in the effective mass (the logarithm of the 
ratio of the correlator on adjacent time-slices) at one particular time-slice for the vr"*", K~^, 
N and S as a function of the number of measurements per configuration. This calculation 
was performed on 664 of the 1194 configurations in the ensemble, those for which we have 
made more than 200 measurements. The correlation functions, after being averaged over the 
sources on each configuration, were blocked in units of 10 configurations (100 trajectories), 
and the uncertainties in the effective mass (EM) on each time-slice were generated with 
the single omission Jackknife procedure. The vr"*" and correlation functions clearly show 
deviations from statistical independence beyond ~ 10 sources per configuration, and by 200 
sources per configuration there is little to be gained by performing additional measurements 
on a configuration. In contrast, measurements of the baryon correlation functions are be- 
having as if they are statistically independent even with 200 sources per configuration. It is 
clear that the statistical uncertainties in the baryon correlators can be further reduced by 
performing even more measurements per configuration. These observations are consistent 
with the arguments regarding the critical source densities. 

An alternative way to investigate this question is to consider the correlation between 
measurements of a correlation function from different sources on the same configuration. A 
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FIG. 3: A log-log plot of the normalized uncertainty in the mean value of the effective mass 
of the TT~^, K~^, N, and H at time-slice t = 34, t = 34, t = 39, and t = 49, respectively, as a 
function of the number of sources on 664 of the gauge-field configurations (those with more than 
200 measurements). The fits correspond to a power-law of the form 6C / (C) = A (Ngy-c)'^. The best 
fit values for the exponent are b = -0.31(2), -0.36(1), -0.51(9), and -0.41(5) for the 7r+, K+, N, 
and H, respectively. (Statistically independent measurements would produce b = -5-) 



natural quantity to consider is an extension of the standard autocorrelator to a source-to- 
source autocorrelator, Xsrc, defined by 



Xsrc(-R; ^o) 



J^C(to;c, s) 



E E E ^(^0' ^2, c)e{su S2 : R) 



- 1 



(3) 



where C{t, c, s) is the correlation function of interest evaluated on time-slice t, configuration 
c and from source s and the function 9{si, S2 : R) is unity if the two sources are separated by 
a 4d-distance |si — S2I < -R. A nonzero value of Xsrc(-R) indicates the presence of significant 
correlations over distances shorter than R. We have calculated Xsrc for a number of the 
correlation functions we analyze but find no sign of deviation from zero even for the case of 
the vr"*". This may in part be due to the poor statistics at small source-source separations 
(see Fig.|2[. 

A further consideration is that for a given number of configurations, at some value of Ns^c, 
the uncertainty in the measurements of a correlation function on a given configuration will 
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FIG. 4: The normalized uncertainties in the measurements of N (left panel) and vr"*" (right panel) 
correlation functions for time-slice t = 10 are shown for some individual configurations as a function 
of the number of measurements on that configuration. The dashed lines significantly below the 
data are the normalized uncertainties on our full ensemble. 

become smaller than the uncertainty in the measurements over the entire ensemble. Once 
this limit is reached, it is pointless to perform further measurements without also increasing 
the ensemble size. Our measurements are far from this limit as is illustrated by fig. |4] where 
the uncertainties in the measurements of tt"^ and N correlation functions on some individual 
configurations are shown as a function of the number of sources and compared to the overall 
uncertainty attained with the full ensemble. 

An important consideration in generating high statistics measurements is the correlation 
between configurations. Ideally, enough trajectories separate each gauge-field in the ensem- 
ble so that they are statistically independent to the precision of the calculation of interest. 
The degree of correlation between configurations dictates the number of measurements that 
should be performed on a given set of configurations before it is more cost effective to enlarge 
the ensemble. In fig. |5]we show the uncertainty at given time-slices of the EM for the tt"*", 
K~^, N and S as a function of the number of gauge-fields on which 50 measurements are 
performed. The configurations are maximally separated in Monte-Carlo time, but an in- 
creasing number of configurations means a reduced separation in Monte-Carlo time between 
each configuration. The curves in fig. |5] correspond to what is expected for statistically 
independent configurations. The 100 maximally separated configurations are separated by 
80 trajectories, the configurations separated by 20 trajectories appear to be contributing 
as one expects for statistically independent configurations (assuming that those separated 
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FIG. 5: A log-log plot of the normalized uncertainty in time-slice t = 34, t = 34 t = 29 and 
t = 39 of the EM of the vr"*", K~^, N and H, respectively, as a function of the number of gauge- 
field configurations, each with 50 measurements. The fits correspond to a power-law of the form 
SC/{C) = A (iVcfg)^ The best fit values are b = -0.52(1), -0.47(1), -0.45(2), and -0.45(1) for the 
vr^ K'^, N, and H, respectively. (Statistically independent measurements would produce b = — 

by 80 trajectories are statistically independent). This is consistent with the hadronic auto- 
correlation times measured on sets of configurations similar to this ensemble in Ref. [TT] . 
f ^ ~ ~ 40. 

D. Computational Costs 

These calculations required significant computing resources to perform; the total cost 
of the measurements was approximately 7 x 10^ JLab 6n cluster node hours (this is an 
older machine with a dual core 3 GHz Pentium D processor per node) distributed over 
various computational facilities. To put this into context, the generation of the gauge-field 
configurations required approximately one-third of this time [H] . Our calculational method 
makes use of hadronic building blocks (partly contracted sets of propagators) which are 
extremely useful for contracting multi-particle correlation functions but inefficient for single 
hadron correlation functions; only 4 x 10^ JLab 6n cluster node hours were directly relevant 
to the calculations presented herein. Nevertheless, it seems that in order to achieve the level 
of precision of the B = 1 correlation functions presented here, propagator generation rather 
than gauge-field generation is the most computationally intensive component of the LQCD 
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calculation. However, even this will be superseded by the calculation of the contractions 
that are required for systems involving more than two baryons (the subject of future work). 



III. EXPECTED SPECTRA 

The form of the correlation functions that are expected to emerge from these calculations 
is a textbook discussion, but is now becoming more relevant as advances in the field are 
enabling more complicated processes to be explored, such as scattering, excited states and 
multi-hadron interactions. Additionally, the accurate statistical sampling we perform in 
the current work brings to light features that have been safely neglected in the past. A 
discussion of the impact of the boundary conditions (here we use anti-periodic temporal 
BCs for the quark fields) on multiple meson correlation functions that were used in a recent 
calculation of K^K^ scattering can be found in Ref. [5], and a more detailed derivation for 
a two particle system can be found in Ref. [IT] . 

For interpolating functions Oa,b, the correlation function that is calculated with anti- 
periodic BCs on the quark-fields is 

Gait) = iTr [e"^^ d\it) 0^(0) 

= 01 0\{^) \k){k\ dsiO) \j) , (4) 

where T is the length of the time-direction and Z = Ti e~^^ is the partition function.^ 

As an example, consider the interpolating operator with baryon number zero, strangeness 
zero (5* = 0), and isospin equal to two (/ = 2) that couples to the tt+tt """-state. This 
state can be written in terms of hadronic field operators as O{0) = ^7r+7r+ vr"'"7r+ + 
Z7r+7r+7r07r0 TT+TT+Tr'^Tr'^ + where the ellipses denote all other possible hadronic field op- 
erators with the same quantum numbers and the Z's are unknown overlap factors. In 
Eq. this operator thus gives non-zero values for (7r~7r~| O{0) |0), (7r~| O{0) |vr+), 
(0| O{0) |7r"'"7r"''), plus all other states with the same quantum numbers as the TT~^Tr~^ 
source. Consequently the corresponding correlation function contains exponentials e~^^ * 
with energies M = Et,+t^+, M = Et,+ — = 0, M = — M = £^^+^,+^,0^,0, 
M = —Et^+t,+t,Ot,o, M = Et^+t,+k+k- ■ ■ ■ In the zero temperature limit, only those ex- 
ponentials with M > Et^+t^+ survive. States with energies less than Et^+t^+ are thermal 
excitations, for instance arising from the process shown in fig. |6| and are quite apparent 
in the measured / = 2 tttt correlation functions and have also been observed in hadronic 
systems involving a static quark [18]. 

Baryon correlation functions are somewhat different, as the interpolating operator for the 
single nucleon, for instance, can couple not only to the N, but also to a state containing the 
N and an even number of vr's, to a p-wave KK state, to a p-wave IIK state, to a p-wave N% 
state, and to any other state with the same quantum numbers as the nucleon. Further, it can 
also couple to backward propagating negative-parity states, such as an s-wave Ntt. Finally, 
the single nucleon interpolating operator can couple forward and backward propagating 

^ Typically, Oa and Ob are closely related; in our calculations, they differ only in the type of smearing of 
the quark fields and in the momentum injection. 
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FIG. 6: A depiction of the thermal contribution to vrvr correlation function. The vertical lines 
indicate the anti-periodic temporal boundaries of the configuration and the grey regions represent 
the 7r^7r+ source and sink. The solid lines correspond to valence quark propagators. 




FIG. 7: A depiction of the thermal A^vr system produced by the single-nucleon interpolating field. 
The vertical lines indicate the anti-periodic temporal boundaries of the configuration and the grey 
regions represent the single nucleon source and sink. The solid lines correspond to valence quark 
propagators, while the dashed lines correspond to a sea quark loop from the gauge-field. 

hadronic states (these are thermal states as they exist only because of the finite temporal 
extent (temperature) of the configuration), an example being a forward propagating N and 
a backward propagating vr or vice versa. These states are simply illustrated by an example 
shown in fig. [7| a iVvr thermal state. Here the finite temporal extent of the configuration 
is indicated by the vertical lines (these should be (anti-) identified). The two grey regions 
correspond to the source and sink interpolating field. In the case depicted, the interpolating 
field at the source is = (uC'j^d'^)!! and that at the sink is A^ = {u^C'^^d)u, suppressing 
spin and color indices. For the usual zero temperature ground state, the source produces 
three valence quarks and the sink annihilates three valence quarks. In the thermal state 
depicted, the source (right grey region) produces two valence anti-quarks and a valence 
quark (solid lines) while also producing, via gluonic interactions, a sea quark-antiquark pair 
(dashed line). The three anti-quarks between the grey regions combine to form an anti- 
nucleon propagating as exp(— M7v(T — t)) where t is the separation between the source and 
sink and we ignore excited states for simplicity. The quark-anti-quark pair propagating 
around the temporal boundary (since two quarks propagate, the boundary appears periodic 
at the hadronic level) contribute a factor of exp(— M^t) where T is the temporal extent of 
the configuration. The resulting contribution to the two point correlator is then 

G{t) ~ Z^^ e-^'-^ e-{M^-M^)t ^ (5) 

corresponding to a state with energy Mn — in the observed spectrum. 

In the case of the correlation function resulting from a single nucleon interpolating opera- 
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tor in the Ai representation of the cubic group, one expects to see a state with energy equal 
to the N mass, Mjv, and also states (a subset of all states) with energy — 2Mt, — 5£'^,r, 
Mat + 6Ej^Tr, and M^r + 2M^ + SEn^^t,, where SE^.^^ is the interaction energy of two tt's in an 
1 = state, SEj^T, is the interaction energy of the Nn system, while SE^^^^ is the interaction 
energy of the A^vrvr system. Particularly disturbing is the state with energy + 6E^tj. 
corresponding to A^vr moving forward in time and a vr moving backwards in time, that con- 
spire to produce a state with an energy that differs from the nucleon mass only by the Ntt 
interaction energy. Such states will be exponentially suppressed by the temporal extent of 
the configuration, however, accurately disentangling such states from the zero temperature 
ground state will ultimately require calculations on ensembles of gauge-field configurations 
with different temporal extents. 

It is important to realize that thermal states are not simply a curiosity that can be 
safely ignored. As we shall see in Section VIII , they dominate the statistical uncertainty of 
baryon correlation functions at large times, providing deviations from the naive form of the 
signal-to-noise ratio. The amplitudes of these states are exponentially suppressed by the 
temporal extent of the configuration times the mass of the backward going hadronic state. 
Consequently, the most important thermal states involve backward propagating pions, and, 
to suppress these states, the product Mj^T must be large. As the chiral limit is approached 
this will become more and more difficult since, in the limit, it is impossible to separate any 
particular state from itself and any number of pions. 



IV. ANALYSIS METHODS 
A. Multi-Exponential Fits 

The high statistics accumulated for this work allows us to perform stable multi- 
exponential fits using a standard minimization. In this section, we explore the deter- 
mination of the ground and excited states as a function of several variables; the number of 
exponentials used in the fit function, Ai'exp, the range of the fit, R, the number of sources 
per configuration, the number of configurations, the blocking time Tbiock, and the (effec- 
tive) anisotropy, ^eff- We present details of our fits for the S, using a correlated fit to the 
smeared-smeared and smeared-point correlation functions. 

To begin, we performed combined multi-exponential fits by minimizing 

x'= Yl iy^(t) - Cs{t)] {cov'X' [yAt') - c:{t')] (6) 

t,t',s,s' 

where i/sit) are the lattice measured correlation functions, s = [SS,SP], and Cov is the 
covariance matrix between both time-slices and correlation functions. The fitting functions 
used are, 

Cssit) = J2 Z'nZ'n e-^"* , Csp{t) = J2 Z'nZ^ e"^"* , (7) 

n n 

where Css (Csp) denotes the smeared-smeared (smeared-point) correlation function. To 
perform these fits we start with a single exponential and perform the correlated fit to Zq , 
Zq and Eq. A selected set of best fit parameters from this fit are used as initial estimates 
for the two exponential fits. This is performed recursively by taking the best fit results from 
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the N exponential fit as an initial estimates to the N+1 exponential fit. With this strategy, 
successful minimizations with up to six exponentials have been performed. However, with 
the inclusion of the fifth and higher exponentials, the minimizer performs poorly, and often 
returns two masses that are degenerate within their uncertainties. Furthermore, as discussed 
in detail in the previous Section, the expected spectrum of states on these anisotropic config- 
urations is such that the resulting masses for the excited states are likely averages of nearby 
energy levels, see also Section IV E below for demonstrations of this. For these reasons, we 
are only confident in the ground state energies extracted in these fits. However, the number 
of exponentials used in a successful minimization plays an important role in minimizing the 
fitting systematic uncertainty. 

The extracted mass of the H as a function of the number of exponentials in the fit form is 
detailed in Table |Tj With the high statistics in this study, fitting a single exponential yields a 
statistical uncertainty of less than 0.2%, with a slightly larger fitting systematic uncertainty, 
however, 50 time-shces must be discarded because of excited state contamination. For 
our multi-exponential fits, the fitting systematic uncertainty is defined to be the standard 
deviation of all successful fits in a given minimization. To define a successful fit, we take a 
fixed length in time, R = tmax — ^min, and a fixed set of initial parameters, and keep all fits 
with an integrated probability distribution Q > 0.1 while varying tmin-^ One observes that 
the statistical and systematic uncertainties are not further reduced by including more than 
three exponentials in the fit. The resulting ground state mass of the S as a function of the 
tmin used in the fit is shown in upper panel of fig. [s] (in a style similar to an effective mass 
plot) with the color and symbol shapes indicating the number of exponentials in the fit. The 
extraction of the nucleon mass is also shown in the lower panel. Increasing the number of 
exponentials in the fit, A'^exp, allows the tmm-interval over which the ground-state energy is 
seen to plateau to be brought closer to the source where statistical uncertainties are much 
reduced. 

The full set of measurements have been used to generate the fits presented in Table |T[ 
and the correlation functions from configurations nearby in Monte-Carlo time have not been 



TABLE I: Multi-exponential fits to smeared-smeared and smeared-point H-correlation functions 
as a function of the number of exponentials, A'cxp- The number of successful fits, Ng^c has been 
defined to be fits with Q > 0.1 for a fixed time-length R with a fixed set of initial parameter values. 
The time window listed is taken as a representative example of a good fit. The first uncertainty is 
statistical while the second is taken from the standard deviation of successful fits, as defined above. 
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^ The quality of fit value, Q, is defined as the integrated probability distribution of with d degrees of 
freedom, Q = dx^'P{x^ ,d), where V{x'^,d) = N{x^)'^^^^^ exp{—x^ /2), with N the normalization 

constant. The lower limit of the integration, Xmim the of the the fit under consideration. 
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FIG. 8: The mass of the H (upper panel) and N (lower panel) extracted from multi-exponential 
fits as a function of the tmin used in the fit. 
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blocked. Blocking is known to be important, since for correlated configurations, unblocked 
correlation functions can lead to underestimates of the true uncertainty. For hadronic quan- 
tities, we expect that the ensemble we have used has an auto-correlation time of about 
40 Monte-Carlo time-steps |TT]. Our calculations have been performed on configurations 
separated by only r = 10. Several different fits were performed to determine the effects 
of blocking on our multi-exponential fits, the results of which are collected in Table [Tlj To 
normalize these fits, a fitting interval, with a range of R = 30, was determined from the 
-^exp = 3 exponential fits. For these fits, the blocking time rbiock has no detectable impact 
on either the statistical or systematic uncertainties (rbiock = 1 corresponds to no blocking, 
while Tbiock = 10 corresponds to blocking every 10 configurations). 



TABLE II: Effects of blocking on the determination of the ground state H mass. 
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The range of time used in the fits also plays an important role in minimizing the uncer- 
tainty. The resulting fits for short and long time-ranges used in three and four exponential 



fits are shown in Table III While the range does not have a significant impact on the statis- 
tical uncertainty, it does significantly reduce the systematic uncertainty in the fit. To have 
such long ranges of statistically useful time-slices, the anisotropy ^ = bs/h, which is 3.5 for 
this ensemble, is crucial. We have not performed calculations with a different anisotropy 
(including isotropy), but this can be qualitatively studied by constructing correlation func- 
tions using only every second or every third time slice, with an effective anisotropy of ^efi = 



1.75 and 1.17, respectively. In the lower section of Table |III[ we display fits of 1, 2 and 
3 exponentials to these reduced sets of measurements. This reduced anisotropy has a sig- 
nificant impact on the resulting uncertainties, particularly for ^cs = 1-17. We were unable 
to find successful four exponential fits, and the number of successful fits with 1, 2 and 3 
exponentials has been reduced. Furthermore, the smaller number of time-slices in the same 
physical extent, reduces our ability to control the systematics of the fits. 

Finally, the impact of the number of sources, A^^src and number of configurations A^^cfg, 
on the uncertainties in the extracted mass of the H has been explored, the results of which 



are collected in Table IV With Ng^c = 50 or A^^fg = 597, the statistical uncertainties with 
-^exp = 3 exponential fits are the same as with the full set of measurements.^ However, in 
both cases the systematic uncertainty is larger than that of the full set of measurements. The 
corresponding dependence of more complicated multi-particle observables on the number of 
configurations and sources are under investigation. 

For multi-exponential fits, it appears that the most important feature in controlling the 
uncertainty the ground state is the number of exponentials with which a successful min- 

^ The set of measurements with varying numbers of sources has been constructed by including all configu- 
rations which have at least Ngic = 100. 
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TABLE III: Effects of fit range, R and anisotropy, ^eff on the determination of ground state E 
mass. 
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imization can be performed. Neither the statistical nor systematic uncertainties improve 
beyond the inclusion of three exponentials in the fits. To have confidence in the iVexp = 3 
exponential fits, the anisotropy is found to be essential. A quantitative exploration of the 
effects of the anisotropy on the stability of multi-exponential fits is desirable, but this would 
be a very costly numerical endeavor. With three or more exponentials, the fitting range, 



TABLE rV: Effects of the number of sources, Nsrc and number of configurations, N^fg on the 
determination of the ground state S mass. 
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number of sources and number of configurations have essentially the expected effect on the 
statistical (and systematic) uncertainties. 



B. Generalized Effective Mass Plots 

Correlation functions on an ensemble of configurations of infinite extent in the time-direction 
become dominated by a single exponential at large times with an argument that is the energy 
of the ground state of the system, 



It is conventional to define the effective mass (EM) from the logarithm of the ratio of the 
correlation function on adjacent time-slices. It is also possible ^ to form a more general EM 
from time-slices separated by t j > 1 



For exponentially decreasing signals with time-independent noise, this will naturally reduce 
the statistical uncertainty in the EM and improve the extraction of energy-eigenvalues as it 
increases the "lever-arm" of the exponential. In such a case, the uncertainty in Mes{tj) 
in Eq. ^ will decrease as 1/tj. Simple correlation functions involving tt's have time- 
independent uncertainties, but this is not the case for baryonic correlation functions, whose 
relative uncertainties grow exponentially with time. We explore the improvements to baryon 
EMs, and ultimately the extraction of baryon masses and the energy-eigenvalues in the vol- 
ume, that result from tj > 1. In fitting an energy to an EM (and other generalizations), 
either the Bootstrap or Jackknife procedures are used to generate the covariance matrix 
associated with the time-slices in the range of the fit.^ This covariance matrix is then used 
to form the x^/dof, which is minimized to determine the energy, and then explored to de- 
termine the uncertainty in this energy. The statistical uncertainty is obtained by finding 
values of the fit parameters where the function attains a value of Xmm + ^■ 

To demonstrate the impact of tj > 1 for baryon EMs, we examine the smeared-smeared 
correlation function of the S-baryon. Figure |9] shows the EMs obtained with tj = 1 (left 
panel) and tj = 10 (right panel). The scatter of the effective mass from time-slice to 
time-slice is significantly reduced with tj = 10 compared with tj = 1, allowing for a clear 
identification of the time range over which it is reasonable to extract the (ground-state) 
mass of the S. Therefore, the systematic uncertainty associated with the fitting range in 
the EM is reduced. The statistical uncertainty in the mass of the S extracted from the EMs 
with the two different values of tj, when fit over the time time-slice interval, are however 



^ This was suggested by K. Jugc in a talk at Lattice 2008, see Ref. [T^, but may have been used eariier. 

^ In the Bootstrap method, iVboot — ^cfg randomly generated bootstrap samples are used after blocking 
over sets of five configurations, while the Jackknife ensembles are constructed by single omission after 
blocking over 10 configurations. We have found consistent results using both methods and by using 
different blockings and values of A'^boot- 



oo 




(8) 



n=0 
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FIG. 9: The left panel shows the conventional EM (i j = 1) from the smeared-smeared H correlation 
function, and the right panel shows the EM for the same correlation function with tj = 10. 



very similar, as can be seen in the resulting fits to time-slices t = 48 to t = 58, 



0.24087 ± 0.00057 ± 0.00080 
0.24060 ± 0.00061 ± 0.00060 



0.68, 
0.44. 



(10) 



The first uncertainty corresponds to the statistical uncertainty in the mass determined from 
the x^/dof minimization, while the second corresponds to systematic uncertainty associated 
with the fitting interval. The systematic uncertainty of this fit is determined by varying the 
fitting interval at each end by 0, ±1, ±2 time-slices, performing a x^/dof minimization over 
each interval and taking half of the spread of the extracted masses. Alternative procedures 
such as using fits to rolling windows of time-slices within the fitting interval return similar 
uncertainties. 

It is interesting to explore how different values of tj modify the form of the covariance 
matrix that is input into the x^/dof minimization. The covariance matrices associated with 
the time-interval t = 48 to t = 51 from these two EMs are shown in Eq. (11 ). They are quite 
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different, with the distant off-diagonal elements becoming more significant for increasing tj. 
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In this comparison, it is important to note that the two extractions make use of different 
parts of the correlation function. The tj = 1 fit uses five time-slices, while the tj = 10 fit 
uses eight well-separated time-shces. 



0.250 rr 
0.248 
0.246 
^ 0.244 - 

■4— ' 
4-1 

0.242 - 
0.240 - 
0.238 - 
0.236 




10 20 30 40 50 60 

t/b, 



0.250 
0.248 
0.246 
3 0.244 

O 

0.242 

-a 

0.240 
0.238 
0.236 



X 



10 20 30 40 50 60 

t/br 



FIG. 10: The left panel shows the conventional EM (tj = 1) from the smeared-point H correlation 
function, and the right panel shows the EM for the same correlation function with tj = 10. 

The EMs from the smeared-point S correlation functions with tj = 1 and tj = 10 
are shown in fig. 10 The scatter in the EM for tj = 1 is substantially less than for the 



smeared-smeared correlation function, as the overlap of the interpolating operator onto the 
ground-state is larger. However, the overlap onto excited states is even larger, and the 
ground state component of the correlation function does not become dominant until later 
times, increasing the fitting systematic uncertainty in the extraction of the ground-state 
mass. 
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C. Prony's Method/Linear Prediction 



As the signal-to-noise ratio of baryon correlation functions degrades exponentially with time, 
it is important to extract the ground-state signal (or excited state signal if that is the state 
of interest) from a range of time-slices starting at the earliest possible time. Significant effort 
has been placed into determining interpolating operators that maximize the overlap onto the 
ground-states of the baryons in order to facilitate this. Further, there has been significant 
effort put into using the variational method [20^ |2T| . for which the correlation functions 
resulting from a number of hadronic interpolating operators are diagonalized on each time- 
slice to give the eigen-energies with the appropriate quantum numbers. A few years ago, 
Fleming suggested that generalizing the EM method to two or more exponential functions 
might be useful in LQCD analysis based on findings of NMR spectroscopists [22] ^- At that 
time, we explored this technique with sets of correlation functions that were available to us 
at that time, and found the method was quite unstable to the statistical fluctuations in those 
measurements. More recently, Lin and Cohen [25] compared this method favorably to the 
variational approach. Given the small statistical uncertainties in the correlation functions 
we are presently considering, and the reduction in the systematic uncertainties achieved with 
tj > 1, we return to explore this technique. 

In LQCD, two-point correlation functions have the form 

Git) = Aoe-"«* + Aie-"i* + ... + Ak-ie-'''-'' + . . . , (12) 

where t denotes the time-slice (time-slices are implicitly taken to be evenly-spaced). It 



follows from Eq. (12) that 



Git + nk) + Gk-iG{t + n{k-l)) + Ck-2 G{t + n{k - 2)) + .. + Co G{t) =0,(13) 

where the integer n is the generalization of tj to the case of a multi-exponential function. 
In order to determine the k coefficients Gi, k equations are required to be formed from the 
measured correlation function. Given the Gi, the roots of 

(e— j'^ + Cfc.i (e^"")'"' + Cfc_2(e-"")'~' + .. + Co = , (14) 

and in particular the a's, provide the energies of the states contributing to the correlation 
function. 



1. One Exponential : The Standard Effective Mass 

In the case of = 1, where the correlation function is assumed to be a single exponential, 
and taking n = 1, 

Git + 1) + Co Git) = , (e-") + Co = , (15) 
and the usual expression for the EM follows trivially. 



^ The method is more generally referred to as Prony's method [23' after Gaspard Riche de Prony who first 
constructed it in 1795 [24 . These techniques and other related methods are known as linear prediction 
theory in the signal analysis community. 
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2. Two Exponentials 



In the case of two exponentials in the correlation function, the most general pair of 
equations that can be used to extract the two effective masses is 

G{t + 2n) + Ci G{t + n) + Co G{t) = 
G{t + 2n + qi) + GiG{t + n + qi) + CoG(t + gi) = , (16) 

where qi is an arbitrary integer off-set between the two equations. Inserting the values of the 
calculated correlation function allows for an extraction of Co,i on each time-slice, j. These 
coefficients are then inserted into 

(e-"")' + Gi (e-"°) + Co = , (17) 

to recover a numerical value for e""°. By choosing n = m = 1, the expressions of Fleming [22] 
are recovered. In order to optimize the two-exponential extraction, a search over values of 
the pair (n, qi) must be performed. A further systematic uncertainty can be assigned from 
this choice. 

The ground-state extracted from the smeared-smeared S correlation function with n = 



gi = 5 is shown in fig. 11 It is clear that the ground-state signal can be isolated from 
the correlation function for a large number of time-slices, many more than using the single 
exponential EM (fig. [9]) alone. We have shown the fit to the ground-state result between 
time slices t = 30 and t = 60. The lower-limit of the time interval was chosen to be within 
an interval for which x^/dof < 1. Extending the fit interval to lower time-slices gradually 



increases the x^/dof, as shown in the right panel of fig. 11, indicating contamination from 
higher energy states. The upper-limit of the fitting interval was chosen to be in the region 
for which backward propagating states (due to the anti-periodic BC's in the time-direction) 
were not visible in the EM (or in the x^/dof). The ground-state S mass we extract from 
this 2-exponential analysis is 

Ms = 0.24109 ± 0.00043 ± 0.00057 , xV^of = 0.38 , (18) 

where the first uncertainty is statistical and the second is the fitting systematic (as de- 
fined previously). The statistical uncertainty in the 2-exponential extraction is significantly 



smaller than that obtained from the one-exponential analysis (Eq. (10)). This is due to the 
substantially increased number of time-slices in the ground-state plateau in the generalized 
EM. 

One aspect of this method that is less appealing is the ambiguity in the association 



of the two roots that result from Eq. (17) to the two states on different time-slices and 



on different jackknife/bootstrap ensembles. This (mis-)identification issue is the cause of 



the anomalously large uncertainties at time-slices 31,. . . ,36 in fig. 11 - this should not be 
interpreted as variance of the signal for the ground state. Additionally, on different time- 
slices and jackknife/bootstrap ensembles, this method can, and likely will, select different 



terms in Eq. (12) particularly for the sub-dominant excited state, adding additional artificial 
variance to the signals for particular energy eigenstates. Consequently, the extracted second 
state is not physically meaningful. 
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FIG. 11: The left panel shows the ground-state extracted from the smeared-smeared H correlation 
function with a 2-exponential Prony determination with n = qi =5, and the correlated fit to the 
time-slices between t = 30 and t = 60. The inner (darker) region corresponds to the statistical 
uncertainty, while the outer (lighter) region corresponds to the statistical and fitting systematic 
uncertainties combined in quadrature. The right panel is the x^/dof for fits between time-slices 
t = tmin and t = 60. 



3. Three and More Exponentials 

The generalization of the method to arbitrary numbers of exponential functions is 
straightforward. In the case of three exponentials, inserting the values of the calculated 
correlation functions, 

G{t + 3n) + C2 G{t + 2n) + Ci G{t + n) + Co G{t) = 
G{t + 3n + qi) + C2G{t + 2n + q,) + CiG{t + n + qi) + Cq G{t + qi) = 
G{t + 3n + q2) + C2G{t + 2n + q2) + Ci G(t + n + ^2) + Co G{t + q2) = ,(19) 

with 5'i 7^ 5'2 7^ allows for an extraction of Co. 1,2 on each time-slice, t. Again, these 
coefficients can be extracted uniquely in terms of the G{t) due to the fact that the system 
is linear. These coefficients Co,i,2 are then inserted into 

(e-"'^)' + C2 (e-"")' + Ci (e-"") + Co = , (20) 
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to recover a numerical value of e~"". Analysis of a given correlation function involves 
searching for the values of the triplet (n, gi, ^2) that optimizes the extraction (in each equality 
in Eq. (19), different n can be used). In this case, statistical fluctuations occasionally result 
in complex roots of Eq. (20) on a particular Jackknife or Bootstrap ensemble. At present, 
we simply omit these contributions in our analysis. Such complex roots correspond to 
an oscillatory solution and arise from short distance noise in the correlation function (or 
nearly degenerate states in the spectrum), and are a well-known issue with the simple Prony 
method. More advanced methods p6] can mitigate this issue, but do not result in improved 
extractions of the ground-state so we do not discuss them in detail. 

The ground-state energy extracted from the smeared-smeared S correlation function with 
n = 10, gi = 3, g2 = 6 is shown in fig. [121 It is clear that the ground-state signal is extractable 
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FIG. 12: The left panel shows the ground-state extracted from the smeared-smeared H correlation 
function with a 3-exponential Prony determination with n = 10, gi = 3,^2 = 6, and the correlated 
fit to the time-slices between t = 10 and t = 52. The inner (yellow) region corresponds to 
the statistical uncertainty, while the outer (red) region corresponds to the statistical and fitting 
systematic uncertainties combined in quadrature. The right panel is the x^/dof for fits between 
the time-slices t = tmm (the horizontal axis) and t = 52. 



from time-slices even closer to the source than with two-exponential analysis. In fig. 12, the 
fit to the ground-state between time slices t = 10 and t = 52 is shown. The extracted mass 
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is 

Ms = 0.24124 ± 0.00032 ± 0.00034 , xV^of = 0.22 , (21) 

with the statistical uncertainty being shghtly less than in the two-exponential analysis. It 
is important to realize that this level of precision corresponds to a statistical uncertainty of 
~ 2 MeV in the S mass. 

We have successfully applied the four- and five- state Prony method to our data but no 
improvement is seen beyond the three-exponential extractions. 



D. Multi-Correlation Function Prony Method 

There are a number of extensions of the Prony method that exist in the literature (see for 
example [2S]), some of which we have investigated in detail. For the correlation functions we 
have in hand, these extensions do not significantly improve on the standard Prony method. 
Typically, these methods are applied in cases where only a single set of measurements is 
available. However, we have two sets of correlation functions (smeared- smeared and smeared- 
point) whose energy spectra are identical in the limit of a large number of configurations. 
It is straightforward to generalize Prony's method to include both correlation functions- 
the matrix-Prony method. This form leads to a further reduction in the uncertainty of 
the extraction of the energy eigenvalues. A similar approach, has been briefly discussed 
in Ref. |27j. 

Assume we have (A^ = 2 in our case) correlation functions from which we want to 
extract the energy levels. If these correlation functions are a sum of exponentials they 
satisfy the following recursion relation, 

My{T + tj)-Vy{T) = ^ , (22) 

where M and V are N x N matrices and and y{t) is a column vector of A^ components 



corresponding to the A^ correlation functions. Eq. (22) implies then the correlation functions 
are 



N 



y{t) = Y,Cnqn\i , (23) 

n=l 

where qn and = exp(m„tj) the eigenvectors and eigenvalues of the following generalized 
eigenvalue problem 

Mq = XVq . (24) 
Given the A^ sets of correlation functions, the masses can be found by determining the 



matrices M and V that are needed in order for the signal to satisfy Eq. (22). Solving 



Eq. (24), then leads to the eigenvalues A„ = exp(m„tj) and the eigenvectors g„ needed to 
reconstruct the amplitudes with which each exponential enters the correlation functions. A 
simple solution can be constructed as follows. First note that 



i^tw t-\-t\Y 



MY,y{T + tj)y{Tf -VY,y{r)y{TY = Q . (25) 
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Clearly, a solution for M and V is 



M 



t+tw 



y{r + tj)y{Tf 



V 



t+tw 

E 

T=t 



-1 



(26) 



where these inverses exist provided that the range, t^, is large enough to make the matrices 
in the brackets full rank {tw > N—1). In our case with two exponentials the range has to be 
two for achieving full rank. Once the eigenvalues, A„ and eigenvectors g„ are determined, the 
amplitudes, C„, can be reconstructed using t as a normalization point. The shift parameter 
tj can be used it improve stability. The above solution is equivalent to determining M and 
V by requiring that 



t+ti 



vl/2 = ^ [Myir + tj) - VyiT)f [My(r + tj) - Vyir)] 



(27) 



is minimized. 

To go beyond extracting two states, one can construct and solve a second order recursion 
relation. The minimization condition of Eq. (27), augmented to contain the second order 
terms in the recursion, can be used to determine the unknown matrices. The resulting 
eigenvalue problem is a second order nonlinear generalized eigenvalue problem which is 
straightforward to solve. However, to isolate the ground-state, which is our present focus, 
the two state model is sufficient and we do not pursue this further. 

To demonstrate how this method works, we return to the S correlation function discussed 



above. Figure 13 shows the generalized EMP for the S mass as a function of time determined 
with a = 2 matrix-Prony extraction, using both the smeared-smeared and smeared-point 
correlation functions. The inset shows the second extracted state in addition to the ground 
state. The extracted value of the S mass, determined by fitting in the time interval t = 11 
to t = 50, is 



0.24097 ±0.00025 ±0.00003 , xV^of = 0.81 



(28) 



The EM of the dominant state in fig. [13] plateaus around time-slice t = 10, and is well- 
defined over a large interval. In addition to being somewhat more visually appealing than 
the previous Prony analyses of single correlation functions, this method provides the smallest 
uncertainties, particularly for the fitting systematic. 

In our final extractions of baryon masses, our EM analysis will use the matrix-Prony 
method. This method yields ground state energies that are in complete agreement with those 
from the other methods discussed. The generalized EMs from the matrix-Prony method 
are consistently clean, and the quality of fits are uniformly good for the ground state. 
Since they involve only one fit parameter, one can easily assess the quality of the fits. 
The procedure for fitting parameters and determining their statistical uncertainty has been 
described in Section IV B[ Systematic uncertainties are calculated by performing fits over 
rolling windows of time-slices within the quoted overall range and looking at the standard 
deviation of the central values of those fits. This is combined in quadrature with a further 
systematic uncertainty that is generated by sampling a large range of possible values of tj 
and tw and taking the standard deviation of the central values of the resulting fits. The 
generalized EMP for the S extracted with the matrix-Prony method for a variety of values 
of tw and tj can be seen in fig. [14 
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FIG. 13: The generalized EMP for the mass of the H using a Matrix-Prony analysis with tj = 7 
and tw = 11) and the correlated fit to the time-slices between i = 11 and t = 50. The inner (darker) 
region corresponds to the statistical uncertainty, while the outer (lighter) region corresponds to 
the statistical and fitting systematic uncertainties combined in quadrature. The inset shows both 
states extracted with the matrrx-Prony method. 



E. Prony-Histograms 

In extracting the energies of the states through the Prony procedure, a set of roots are 
produced on each time-slice for each member of the Bootstrap or Jackknife ensemble. In 
general these roots are real and there is an ambiguity in associating the roots with energy- 
levels in the finite volume (only the single particle masses are approximately known). In 
order to aid identification of energy-levels it is useful to form histograms of the complete 
set of roots generated through the Bootstrap procedure. The simplest histogram is formed 
by accumulating all of the roots obtained on a sub-set of, or all of, the time-slices over all 
bootstrap/jackknife ensembles. The dominant components of the correlation function will 
appear as well-defined peaks in the histogram. 

In most cases, this histogramming procedure produces very similar results when either 
two, three or four exponential Prony, or matrix-Prony analyses are used. Only atypically do 
the higher exponential analyses reveal a clean state that is not present in the two-exponential 
analysis. Additionally, since our baryon correlation functions are asymmetric in time because 
of the parity projectors used in Eq. (|2]), noise is reduced in these histograms by separately 
accumulating the roots over the two half configuration. As expected, the excited states have 
a larger presence in the smeared-point correlation function. An example histogram is shown 



in fig. 15, corresponding to the S correlation function analyzed in fig. 13 There is one 
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FIG. 14: The generalized EMP for the mass of the S using a Matrix-Prony analysis for a variety 
of values for tj and tvK- 



clear peak in the histogram, corresponding to the S ground-state and one broad structure 
at higher mass, which the histogram suggests is hkely to be a collection of closely spaced 
states that currently are not resolvable. This interpretation of the excited state is consistent 
with expectations for the S spectrum and with the instability of the extractions of the first 
excited state in the exponential fits discussed earher. 



V. MESON SPECTRUM 

The 77"*" correlation functions do not suffer from exponential signal-to-noise degradation for 



configurations of infinite temporal extent (in Section VIII , we will find that this is not true 



for a finite time-direction even for the pion). As a result, they can be calculated with small 



statistical uncertainty on each time-slice, as shown in fig. 16 As the EM for the vr"*" does 
not exhibit a plateau, the 7r+ mass is determined by fitting cosh (M7r(t — -|)) to a (large) 
number of time-slices of the correlation function. Performing a double cosh fit to time-slices 
t = 21 to t = 41 yields 

= 0.06936 ± 0.00012 ± 0.00005 , xV^of = 0.73 , (29) 

where the first uncertainty is statistical and the second is fitting systematic. The statis- 
tical uncertainty in the mass is determined with the Jackknife procedure, and the fitting 
systematic is determined by varying the fitting interval over a reasonable range. 

The second set of peaks that are visible in the histogram are at an energy consistent with 
the 1 = 1 KKn state (with a threshold at -|- 2Mk ~ 0.2636) that can couple to the 
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FIG. 15: The histogram of the positive roots extracted from time-slices t = 11 to t = 50 from 
N = 2 matrix-Prony analysis of the H correlation functions, with tj = 7 and t\\r = 11. 



source that produces a single vr"*". With even greater statistics, the energy of this state could 
be calculated with enough precision to extract the I = \ Ktx and / = KK scattering 
lengths and the 1 = 1 KKtt three-body interaction. An expression for the energy-levels of 
this system in a finite volume in terms of the KK and Kn scattering amplitudes and various 
three-body interactions has recently been derived [28] and would be useful in analyzing this 
state. There is no clear peak that can be associated with 7 = 1 TTTrvr, which one would 
naively thought would have been present. It must be the case that the source does not 
couple with any appreciable strength to this state. 



The EM associated with the smeared-point kaon correlation function is shown in fig. 17 
along with the bootstrap-Prony histogram. Despite the appearance of the EM, no plateau 
is found in the EM, and the kaon mass is extracted by fitting cosh (Mi^(t — |^)) to a number 
of time-slices of the correlation function. Performing a double cosh fit over the time-slices 
between t = 29 and t = 49, yields a K~^ mass of 

Mk = 0.097016 ±0.000099 ±0.000033 , xV^of = 1.01 . (30) 

The excited state(s) that are seen in the histogram in fig. n/nare consistent with the I = \ 
KKK. A better measurement of this state, in analogy witnthe pion correlation function, 
would allow for a determination of the 1 = KK scattering amplitude. 
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FIG. 17: The upper panel shows the EM for the smeared-point correlation function, while 
the lower panel shows the associated Prony histogram. 
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VI. GROUND-STATE BARYON SPECTRUM 



With the methodology we have presented in Section IV, we are in a position to extract 
the masses of the lowest-lying octet baryons. The S correlation functions have been used 
extensively to demonstrate the strengths and weaknesses of the various methods, with the 



resulting mass extraction given in Eq. (28), and we do not repeat that discussion here. 

The matrix-Prony method applied to the smeared-smeared and smeared-point correlation 
functions associated with the E, A and produces the Prony-histograms and generalized 



EMs shown in fig. 18, fig. 19, and fig. 20, Fitting the S EM between time-slices t = 12 to 



t = 47 yields S mass, 

Ms = 0.22811 ±0.00028 ±0.00018 , xV^of = 0.77 . (31) 
Fitting the A EM between time-slices t = 12 to t = 52 yields A mass. 

Ma = 0.22255 ±0.00028 ±0.00005 , xV^of = 1.21 . (32) 
Finally, fitting N the EM between time-slices t = 11 to t = 40 yields N mass, 

Mn = 0.20682 ± 0.00032 ± 0.00010 , xV^of = 1.5 . (33) 

The results of the best extractions of the ground-state baryon masses using multi- 
exponential fitting and the matrix-Prony method, which give consistent results for each 
species of baryon, are collected in Table |Vj These results are completely consistent within 
their uncertainties, giving us confidence that our extractions are correct. 



TABLE V: The ground-state masses of the = 2 baryons extracted by fitting four exponentials 
and by the matrix-Prony method. The first uncertainty is statistical while the second is the fitting 
systematic. 
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VII. NEGATIVE-PARITY EXCITED BARYON STATES 

The interpolating operators that produce even-parity baryons moving forward in time also 
produce negative-parity partners moving backwards in time. As the interpolating operators 
couple to continuum states such as A^vr, it is possible that, by using Liischer's method 
(and ideally, multiple spatial volumes), the phase-shifts for meson-baryon scattering can be 
extracted in channels with contributions from disconnected diagrams. 

In addition to excited single baryon states, and the continuum states that carry zero units 
of momentum in the volume, there are also continuum states where each hadron carries one 
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FIG. 18: The upper panel shows the generalized EM for the mass of the S using a Matrix-Prony 
analysis with tj = 9 and tw = 17, and the correlated fit to the time-slices between t = 12 and 
t = 47. The inner (darker) region corresponds to the statistical uncertainty, while the outer (lighter) 
region corresponds to the statistical and fitting systematic uncertainties combined in quadrature. 
The inset show the same ground-state EM plot along with that of the excited state (light points). 
The lower panel shows the associated Prony histogram of the positive roots for the time-slices 
t = 12 to t = 47. 34 
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FIG. 19: The upper panel shows the generalized EM for the mass of the A using a Matrix-Prony 
analysis with tj = 9 and twr = 11, and the correlated fit to the time-slices between t = 10 and 
t = 47. The inner (darker) region corresponds to the statistical uncertainty, while the outer (lighter) 
region corresponds to the statistical and fitting systematic uncertainties combined in quadrature. 
The inset show the same ground-state EM plot along with that of the excited state (light points). 
The lower panel shows the associated Prony histogram of the positive roots for the time-slices 
t = 10 to t = 47. 35 




FIG. 20: The upper panel shows the generahzed EM for the mass of the N using a Matrix-Prony 
analysis with tj = 7 and tw = 2, and the correlated fit to the time-slices between t = 11 and t = 40. 
The inner (darker) region corresponds to the statistical uncertainty, while the outer (lighter) region 
corresponds to the statistical and fitting systematic uncertainties combined in quadrature. The 
inset show the same ground-state EM plot along with that of the excited state (light points). The 
lower panel shows the associated Prony histogram of the positive roots for the time-slices t = 11 
to t = 40. 36 



or more units of momentum in the volume, while having vanishing total momentum. The 
lowest energy state containing hadrons A and B with back-to-back momenta ±p = ±^n 
(where n is an integer triplet) occurs at 



EZ-^|Ml+[^-^) . (34) 

In attempting to unravel the spectrum of states contributing to the correlation functions, 
we must also consider such continuum states. 

The lowest-lying negative parity state that is expected to couple to the interpolating 
operator for the single nucleon is the s-wave iVvr state (more precisely, we refer to the Af 
representation of the hyper-cubic group), which has a threshold, neglecting interactions, of 



M^ + Mn = 0.27618±0.00034±0.00011. Fitting the EM shown in fig. [21] between time-shces 
t = 93 to t = 119 yields, 

En^ = 0.2861 ± 0.0011 ± 0.0020 , xV^of = 0.91 , (35) 

significantly above threshold. Therefore, we conclude that this state is an s-wave ttN scat- 
tering state with isospin J = |, as it has energy considerably below that of the first mo- 
mentum excitation in the volume at eJ^J^^ = 0.33889 ± 0.00042 ± 0.00014, and the vrvrA^ 
state. Given that there are no other channels that are energetically allowed for this state 
to mix with, the s-wave vrA^ phase-shift^ in this channel can be extracted at an energy of 
6E^N = ^^(1-) -Mn-M^ = 0.0095 ± 0.0011 ± 0.0020 {6E^n = 15.3 ± 1.8 ± 3.2 MeV) 
where the uncertainties are dominated by the uncertainty in E^^i^y It is important to note 

that this channel receives contributions from disconnected diagrams, and in the calculation 
we are doing, these contributions are completely accounted for in the gauge-configurations. 
Using the standard Liischer procedure, a phase-shift of 6t,n = —26 ± 7 ± 6 degrees is found 



at this energy. The Prony histogram in fig. 21 shows significant structure in this channel. 



and one could argue that there is a single level at -E ~ 0.45, but this would require further 
exploration. 

For the negative-parity state that couples to the interpolating operators for the A, the 
situation is not so clean. The thresholds for the lowest-lying continuum states, Svr and 
NK, are located at + = 0.29747 ± 0.00030 ± 0.00019 and Mn + Mk = 0.30384 ± 
0.00033 ± 0.00011, respectively, in the absence of interactions. The corresponding lowest- 
lying states with one unit of momentum occur at E^^l~^ = 0.35857 ±0.00037 ±0.00023, and 
eJ^^^ = 0.35763 ± 0.00039 ± 0.00012. Fitting the EM shown in fig. 
t = 88 to t = 117 yields. 
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between time-slices 



Ej^^r-^ = 0.2983 ±0.0008 ±0.0004 , xV^of = 1.02 . (36) 

This is, within uncertainties, at the threshold for Stt or NK. The eigenstates will be a 
combination of these two systems and it is likely that we have not resolved the two nearby- 



states in the EM, and the result in Eq. (36) is actually an average of two closely-spaced 
energies. 

® Here we ignore possible contributions from L = 4,6, .. . partial waves that also contribute in the A'^ 
representation of the cubic group H(3). 
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FIG. 21: The upper panel shows the generalized EM for the lowest-lying negative-parity state 

coupling to the N-source using a Matrix-Prony analysis with tj = 5 and tyy = 8, and the correlated 

fit to the time-slices between t = 93 and t = 119. The inner (darker) region corresponds to the 

statistical uncertainty, while the outer (lighter) region corresponds to the statistical and fitting 

systematic uncertainties combined in quadrature. The lower panel shows the associated Prony 

histogram of the positive roots for the time-slices t = 93 to t = 119. 
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FIG. 22: The upper panel shows the generahzed EM for the lowest-lying negative-parity state 

coupling to the A-source using a Matrix-Prony analysis with t j = 1 and tw = 5, and the correlated 

fit to the time-slices between t = 88 and t = 117. The inner (darker) region corresponds to the 

statistical uncertainty, while the outer (lighter) region corresponds to the statistical and fitting 

systematic uncertainties combined in quadrature. The lower panel shows the associated Prony 

histogram of the positive roots for the time-slices t = 88 to t = 117. 
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FIG. 23: The upper panel shows the generahzed EM for the lowest-lying negative-parity state 

coupling to the S-source using a Matrix-Pro ny analysis with tj = 7 and tw = 2, and the correlated 

fit to the time-slices between t = 95 and t = 118. The inner (darker) region corresponds to the 

statistical uncertainty, while the outer (lighter) region corresponds to the statistical and fitting 

systematic uncertainties combined in quadrature. The lower panel shows the associated Prony 

histogram of the positive roots for the time-slices t = 95 to t = 118. 

40 



For the lowest- lying negative-parity state (s) produced by the interpolating operator for 
the E, the situation is even more complicated. The threshold of the non- interacting s- 
wave Ett state is at Ms + = 0.29747 ± 0.00030 ± 0.00019, for the Avr state is Ma + 
M^ = 0.29191 ± 0.00030 ± 0.00007, and for the NK state is + Mk = 0.30384 ± 
0.00033 ± 0.00011. Therefore, in this large volume, we expect to observe three eigenstates 
that are nearly degenerate. The lowest-lying states with one unit of momentum occur at 
E^"l=^ = 0.35857 ± 0.00037 ± 0.00023, = 0.35341 ± 0.00030 ± 0.00007, and E^^^^ = 

0.35763 ± 0.00039 ± 0.00012 and are well separated from the n = states. Fitting the EM 



shown in fig. [23] between time-slices t = 95 to t = 118 yields, 

E^^i-^ = 0.3068 ±0.0011 ±0.0011 , xV^of = 0.80 , (37) 

in the region where one expects to find three closely-spaced states, corresponding to the 
eigenstates dominated by Svr, Att, and NK. Given how closely spaced these states are 



expected to be, the extraction in Eq. (37) is likely a complicated average of three energies. 

The situation is no better for the lowest-lying negative-parity states that are expected 
to couple to the interpolating operator for the S. The lowest-lying s-wave continuum states 
are vrS, AK and TiK. The threshold for these states, in the absence of interactions, are 
e["^=° = 0.31033 ± 0.00028 ± 0.00006, sj^f" = 0.31957 ± 0.00030 ± 0.00006, and e\'^^=° = 
0.32513 ± 0.00030 ± 0.00018, respectively. The corresponding states where both hadrons 
carry one unit of momentum have thresholds state = 0.37058 ± 0.00033 ± 0.00007, 

e["^=i = 0.37214 ± 0.00035 ± 0.00007, E^-^j^^ = 0.37731 ± 0.00034 ± 0.00021, respectively. 
Therefore, we expect to observe two sets of three nearly degenerate eigenstates. Fitting the 



EM shown in fig. 24 between time-slices t = 91to)f: = 118 yields. 



E i_ = 0.3243 ±0.0010 ±0.0009 , xV^of = 0.72 , (38) 

in the region where one expects to find three closely-spaced states, corresponding to the 
eigenstates dominated by n = Hvr, AK, and T,K. Given how closely spaced these states 



are expected to be, the extraction in Eq. (38 ) is likely an average of three unresolved energies. 
There are hints of a couple of other peaks in the Prony-histogram, but nothing conclusive. 

The results of the best extractions of the ground-state baryon masses using multi- 
exponential fitting and the matrix-Prony method, which give consistent results for each 



species of baryon, are collected in Table VI 



TABLE VI: The masses of the lowest-lying = 2 states with unit baryon number extracted 
by fitting three exponentials and by the matrix-Prony method. The first uncertainty is statistical 
while the second is the fitting systematic. 





Exponential Fitting 


Matrix-Prony 


state 


btM 


range 


xVdof 
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btM 


range 


xVdof 




0.2871(18)(10) 


90-117 


1.11 


0.28 


0.2861(11)(20) 


93-119 


0.91 


A(r) 


0.2954(05)(15) 


90-113 


0.89 


0.64 


0.2983(8)(4) 


88-117 


1.02 


s(r) 


0.3074(15)(15) 


90-118 


1.02 


0.41 


0.3068(11)(11) 


95-118 


0.80 




0.3261(09)(15) 


89-115 


1.07 


0.36 


0.3243(10)(9) 


91-118 


0.72 
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FIG. 24: The upper panel shows the generahzed EM for the lowest-lying negative-parity state 

coupling to the H-source using a Matrix-Prony analysis with tj = 5 and tw = 11, and the correlated 

fit to the time-slices between t = 91 and t = 118. The inner (darker) region corresponds to the 

statistical uncertainty, while the outer (lighter) region corresponds to the statistical and fitting 

systematic uncertainties combined in quadrature. The lower panel shows the associated Prony 

histogram of the positive roots for the time-slices t = 91 to t = 118. 
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VIII. SIGNAL-TO-NOISE RATIOS 



Many observables of importance to particle physics that are currently being calculated with 
LQCD, such as the pion decay constant and the Gasser-Leutwyler coefficients, require the 
calculation of mesonic correlation functions. Statistical fluctuations on each time-slice of 
these correlation function are well-behaved. In contrast, as argued by Lepage [10], correlation 
functions involving one or more baryons exhibit exponentially growing statistical noise. In 
the case of a single positive parity nucleon, the correlation function has the form 



{e^{t)) = J2 rf (iv"(x,t)iv^(o,o)) ^ z, 



(39) 



where N is an interpolating field that has non-vanishing overlap with the nucleon and the 
angle brackets indicate statistical averaging over measurements on an ensemble of configu- 
rations. The variance of this correlation function is 



y rj«r7^^^a(^^ ^)^/3(y^ ^)^.^o^ 0)^^(0^ 



x,y 



^Stt e — Zjj^ e — i> Zstt e 
and therefore, as Lepage [10] argued, the noise-to-signal ratio behaves as 

More generally, for a system of A nucleons, the noise-to-signal ratio behaves as 

a 1 



..AM. 



X 



N 



(40) 



(41) 



(42) 



-AM Alt 



the noise-to-signal associated 



Therefore, in addition to the signal itself falling as G ~ e 
with the correlation function grows exponentially, as in Eq. (42). 

These arguments are constructed for a system with an infinite time-direction and are 
modified in an important way for systems with a finite time-direction with given BCs. 
The calculations that are presented in this work have employed anti-periodic BCs in the 
time-direction. With such BCs the positive parity nucleon correlation function in Eq. (39) 
becomes 



{eM{t)) ^ Zn 6-"'^' + Zm. e+^^-(*-^) 



(43) 



where E^-k is the energy of the lowest-lying negative-parity state in the volume, which, for 
this ensemble of configurations, is a continuum nucleon and pion at rest. The arrow denotes 
the behavior of the correlation function far from source (in both time-directions). Further, 
the correlation function dictating the behavior of the variance of the nucleon correlation 
function is modified similarly, with Eq. (40) becoming 
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+ A^ e 
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+ Ao e-*'^^ + 



(44) 
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The first term in Eq. (44) arises from Svr's propagating forward and Stt's propagating back- 



wards, the second term arises from 2'7r's propagating forward along with one vr propagating 
backward and vice versa, the third (time-independent) term arises from a nucleon prop- 
agating forward and a nucleon propagating backward, and the ellipses denotes terms in- 
volving larger masses. As the negative-parity state is more massive than the nucleon, the 



nucleon is the dominate component in the correlation function, Eq. (43), for a number of 
time-slices beyond the mid-point of the configuration. From this argument, one expects 
to see the signal-to-noise ratio degrade even more rapidly than the expectation shown in 
Eq. (41) in the time-slices near the mid-point of the configuration where the correlation 
function is still dominated by the nucleon. One expects to find regions of the correlation 
function, depending on the structure of the source, which have the noise-to-signal scaling as 



and e'^p^^ or combi- 



nations thereof. 

The high statistics calculations we are presenting here enable a detailed study of the 
behavior of the signal-to-noise ratio associated with the correlation functions formed with 
quark propagators generated with anti-periodic BCs. It is useful to form the effective 
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FIG. 25: 
with tj - 

noise-to-signal ratio, as defined in Eq. (46). The horizontal lines correspond to the energy scales 
1 



The left panel shows the EM of the smeared-point N correlation function formed with 
= 3. The right panel shows the energy-scale, Eg, associated with the growth of the 
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FIG. 26: The left panel shows the EM of the smeared-smeared N correlation function formed with 
tj = 5. The right panel shows the energy-scale, Es, associated with the growth of the noise-to- 



signal ratio, as defined in Eq. (46 ). The horizontal lines correspond to the energy scales nip — iM^r, 



^Mt^, nip, nip + ^Mtt, and nip + ^M.^ (from lowest energy to highest energy). 



noise-to-signal plot, in analogy with the EMs. On each time slice, the quantity 



S{t) 



x{t) 



(45) 



is formed, from which the energy governing the exponential behavior can be extracted via 



tj 



log 



Sit + tj) 
S{t) 



(46) 



If the correlation function is dominated by a single state, and a single energy-scale determines 
the behavior of the noise-to-signal ratio, the quantity Es(t; tj) will be independent of both 
t and tj. 



In fig. 25, the full EM of the smeared-point nucleon correlation function is shown (with 
tj = 3), and in fig. 26, the full EM of the smeared-smeared nucleon correlation function 
is shown (with tj = 5). Also shown are the energy-scales associated with the growth of 



the noise-to-signal ratio from Eq. (46), with uncertainties generated using the Jackknife 



procedure. Considering the smeared-point correlation function in fig. 25 after time-slice 
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t = 35 or so, the correlation function is dominated by the ground-state nucleon which 
persists until time-slice t ~ 70. Beyond this time-slice the backward propagating negative- 
parity A^vr-state becomes dominant. Between time-slices t ~ 40 and t ~ 50, the noise- 
to-signal ratio is determined by the expectation of rup — |M^. However, after t ~ 50 the 
signal-to-noise ratio degrades exponentially faster than this, and by t ~ 65 the relevant 
energy-scale is ~ nip + \M.,^ and increasing with t. Similar behavior is clear in the smeared- 
smeared correlation function, for which the nucleon ground state dominates from an earlier 
time-slice. 

It is clear from this analysis of the noise-to-signal ratio, that the length of the time- 
direction of these configurations and resulting thermal states are limiting the precision of 
the ground-state nucleon mass determination. This will be even more true for the multiple 
baryon correlation functions for which the signal-to-noise degrades exponentially faster than 
in the single nucleon correlation functions. Increasing the length of the time direction will 
lead to exponential improvement of the correlation function at large times where the nucleon 
component dominates the correlation function. It is interesting to note that the coefficients of 
the backward propagating contributions to the noise-to-signal ratio are suppressed by powers 
of e-2^'rT_ Qj-^ ^YiQ current configurations with T = 128, in order to reduce the contribution 
to the noise from the nip — ^M^r component by an order of magnitude, the time extent would 
need to be increased to T ~ 192. This will reduce the nip + component by a factor 
of ~ 84 and the nip + iM^r component by ~ 770. Such an increase in the temporal extent 
would significantly decrease the statistical uncertainties with which ground-state signals are 
extracted. 

The noise-to-signal analysis of the S correlation functions is somewhat more complex, 
because there are a number of low-lying states which can contribute to the variance. For 
the SS noise correlation function, the lightest intermediate states that can couple to quark- 
content ssussu are KKt], rjrjTi and 777777. The full EMs and the Eg plots for the smeared- 



point and smeared-smeared S correlation functions are shown in fig. 27 and fig. 28 In both 



the smeared-point and smeared-smeared S correlation functions, the noise-to-signal ratio is 
growing exponentially slower than naive expectations, until about time-slice t ~ 55. As the S 
ground-state dominates the smeared-smeared correlation function beyond t ~ 40, this allows 
for a extraction of the mass with higher precision than expected. This suggests that the noise- 
source does not couple to the low-lying mesonic states as strongly as expected, and that more 
massive mesonic states are dominating the noise over many time-slices. However, eventually, 
for t ^ 55, the growth of noise overshoots the original Lepage expectation (indicated by the 



lowest horizontal line in figs. 27 and 28 



Whilst we are primarily interested in noise in the baryonic sector, it is interesting to 
note that the mesonic correlation functions also suffer from similar issues. According to 
the above arguments, the pion correlation function on lattices at zero temperature (infinite 
temporal extent) will have noise that is independent of time (up to fiuctuations) while the 
kaon will have noise that grows exponentially with the small energy difference niK — \t^-k — 
^niri. However at finite temperature, the noise correlation functions of both systems receive 
additional contributions that grow faster than the above expectations. This is shown in 
fig. [29} 
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FIG. 27: The left panel shows the EM of the smeared-point H correlation function formed with 
tj = 3. The right panel shows the energy-scale, Es, associated with the growth of the noise-to- 
signal ratio, as defined in Eq. (46). The horizontal lines correspond to the energy scales — |m.^, 
1 



ms - Mk 
m~ + m„ - 



rrir 



Mk + 



\Mt^ (from lowest energy to highest energy). 



and 



IX. SCALING WITH COMPUTATIONAL RESOURCES 



An important component of our current work is to address the future requirements for 
LQCD calculations in nuclear physics, a field characterized by small energy scales in heavy 
systems, for example, the 2 MeV binding energy of the ~ 2 GeV deuteron. In fig. 30 , we show 
the extracted mass of the vr"*", , N and S as a function of the number of configurations 
in the ensemble for both the exponential and matrix-Prony analysis methods. The full set 
of measurements performed on each configuration are included, and the fitting intervals are 
chosen to optimize the extraction for each ensemble size. In each case, the uncertainty in the 
mass is reduced, as expected, with increasing ensemble size, and the mass extracted from 
the smaller ensembles tends to be less than that from the larger ensembles. Fig. 31 shows 
the fractional uncertainty in the mass of the 7?+, , N and S, associated with the results 
in fig. |30| as a function of the number of configurations. An extrapolation can be performed 

of the form 5M/M 



with a fit to the uncertainties in fig. 



31 



in these fits are -0.55(4), -0.51(3), -0':^(4), -0.67(6) for the 7r+, K+ , N and S, respectively. 



The exponents extract 
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FIG. 28: The left panel shows the EM of the smeared-smeared H correlation function formed with 
tj = 3. The right panel shows the energy-scale, Es, associated with the growth of the noise-to- 
signal ratio, as defined in Eq. (46). The horizontal lines correspond to the energy scales — |m.^, 
1 



ms - Mk 
m~ + m„ - 



rrir 



Mk + 



\Mt^ (from lowest energy to highest energy). 



ni 



H + ^Mtt, m-s, and 



The dependence of our results for hadron masses on the number of sources used in the 
calculations is explored in fig. 32 where we show the fractional uncertainty in the mass of 
the TT"*", K~^, N and S as a function of the number of sources used on each configuration. In 
this figure we use an ensemble of 1012 configurations, those on which there are at least 100 
measurements. A simple fit of the form 6M/M = AN^^.^ returns exponents b = —0.03(2), 
-0.65(19), -0.41(3), -0.40(6) for the 7r+, K+, N and S, respectively. 

The results of this analysis can be simply summarized for baryons (averaging over the 
nucleon and S) as 

~ ' (47) 

For mesons, a similar scaling is seen, with a somewhat worse scaling with the number of 
measurements per configuration in the case of the pion, consistent with the saturation seen in 
fig. |3| This functional form enables us to quantify the relative benefit of increasing the num- 
ber of sources per configuration compared to increasing the total number of configurations. 
The costs involved in this are as follows: 
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FIG. 29: The energy scale, Eg, associated with the growth of the noise-to-signal ratio in the vr"'" 
(left) and (right) smeared-smeared correlation functions using tj = 3. The horizontal lines 
correspond (from lowest to highest) to 0, Mj^ for the pion and Mk — \'^ri — \'^K for the kaon. 

• Gauge configuration generation: The total cost of generating the ensemble of 1194 
gauge configurations was 2M JLab-6n cluster node-hours and the production took a 
significant amount of wall-clock time. Configuration generation costs scale linearly 
with the number of configurations once a Monte-Carlo trajectory has thermalised (in 
this case the overhead of thermalisation was approximately 10%). In order to generate 
significantly larger ensembles (containing 10^ or 10^ gauge fields) in a reasonable wall- 
clock time, it will be necessary to run multiple trajectories in parallel. Given wall-clock 
time and memory constraints, an individual trajectory will produce (9(1000) gauge- 
field configurations that are useful for measurements. Consequently the thermalization 
overhead will conservatively remain at about 10%. Each configuration requires ~ 
2 X 10^ JLab-6n node-hours to produce. 

• Measurement calculations: The total cost of computing all of the measurements per- 
formed in this work was 7 x 10^ Jlab-6n node-hours. The cost to generate the 245 
light-quark and strange-quark propagators per configuration on the 1194 configura- 
tions in this ensemble was ~ 3M Jlab-6n node-hours, while the cost to generate the 
baryon and meson blocks (used at intermediate stages of the calculations) was ~ 3.5M 
Jlab-6n node-hours. Contracting the blocks to accomplish the desired measurements 
(one, two, ... baryons, one, two,... mesons and so forth) cost ~ 0.5M Jlab-6n node- 
hours. If propagators on a given configuration are computed in sets of 100, the initial 
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FIG. 30: The extracted masses of the 7r+, N and H as a function of the number of configura- 
tions (with the full set of measurements performed on it). Statistical and systematic uncertainties 
have been combined in quadrature. For the baryon states, both the matrix-Prony and exponential 
fits are shown. 



overhead of constructing deflation vectors in the EigCG algorithm becomes negligible 
(at the 1% level) and can be eliminated for further sets of calculations by storing the 
eigenvector information. On typical machines, each set of propagators and associ- 
ated hadron blocks (technically, not an efficient way to calculated the single hadron 
spectrum, but critical for two and more hadron calculations) requires 22 JLab-6n 
node-hours to produce. 

Anisotropy: The anisotropy of the lattices us ed in our calculations proved useful in re- 



ducing systematic errors in our fits (see Table III), providing approximately a 1/ ^/C 



reduction. However, the cost of producing gauge-field configurations and propagators 
scales as approximately for the same physical extent (one power arises from the ad- 
ditional time-slices and one power arises from the worsening condition number of the 
Dirac operator). Comparing these exponents, we would conclude that using anisotropic 
configurations is not ideal. However for more complicated multi-hadron systems where 
useful fit ranges are much reduced in physical units, the anisotropy will likely prove 
to be very useful. This remains to be investigated further in subsequent studies. 



Using this information and the scalings in Eq. (47), we can address the question of how 
much computation is required to achieve a particular level of statistical precision. With the 
current data this is only possible in the single hadron sector; ongoing analyses will address 
the i? > 1 in the near future. To halve the uncertainty in the determinations of the ground- 
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FIG. 31: The fractional uncertainty in the extracted masses of the 7r+, K~^, N and H as a function of 
the number of configurations (with the full set of measurements performed on it) for the exponential 
analysis. Statistical and systematic uncertainties have been combined in quadrature. The curves 
correspond to fits of the form 6M/M = AN^!^^. The exponents extract in these fits are -0.55(4), 
-0.51(3), -0.38(4), -0.67(6) for the 7r+, K+, N and H, respectively. 



state baryon masses (calculating the nucleon mass at the ~ 1 Me V- level), an increase in 
the number of configurations by a factor of four, or of the number of measurements per 
configuration by a factor of 5.6 is required. Achieving this precision by performing more 
measurements on the existing set of configurations (~ 1100 additional measurements on 
1200 configurations) would cost 30 M JLab-6n node-hours. Achieving the same precision by 
generating an additional configurations and performing the same number of measurements on 
them (3600 configurations with 245 measurements) would require 27 M JLab-6n node-hours. 
Both approaches have similar cost at this level of precision, but for further improvements, 
the generation of additional configurations will be more efficient. Additionally, the second 
approach will further improve the uncertainty for observables such as the pion mass that 
have saturated in terms of the number of sources per configuration. This approach would 
clearly be of more benefit to the broader community. 



X. CONCLUSIONS 

The energy-scales that arise in nuclear physics are typically in the MeV range, and in order 
for LQCD to have significant impact in this field, baryon masses (and energy-eigenstates in 
the volumes relevant to scattering processes) must be calculable with uncertainties that are a 
fraction of an MeV (including isospin-breaking and electromagnetic interactions, quark mass, 
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FIG. 32: The fractional uncertainty in the extracted masses of the vr"*', ^ N and H as a func- 
tion of the number of sources used on each configuration (1012 configurations were used in this 
study) for the exponential analysis. Statistical and systematic uncertainties have been combined 
in quadrature. The curves correspond to fits of the form 5M/M = AN^^^. The exponents extract 
in these fits are -0.03(2), -0.65(19), -0.41(3), -0.40(6) for the 7r+, K+ , N and H, respectively. 



lattice volume and lattice spacing extrapolations). Current computational resources do not 
permit such calculations. In this work we have performed the first high-statistics study of 
baryon correlation functions to better understand a number of issues that will impact the 
precision with which quantities of importance to nuclear physics can be determined with 
LQCD. In the future, we will extend our analysis to look at observables in the B > 1 baryon 
sector. 

At the single lattice spacing, lattice volume and unphysical light quark mass used in this 
work, we find the following set of ground state masses 





= 390.3(0.7)(0.3)(2.5) MeV, 


Mk 


= 546.0(0.6)(0.2)(3.6) MeV, 


Mn 


= 1163.9(1.8)(0.6)(7.6) MeV, 


Ma 


= 1252.4(1.6)(0.3)(8.2) MeV, 


Ms 


= 1283.7(1.6)(1.0)(8.4) MeV, 


Ms 


= 1356.1(1.4)(0.2)(8.8) MeV, 


En(1/2-) 


= 1610(06)(11)(11) MeV, 


-Sa(1/2-) 


= 1679(05)(02)(11) MeV, 


-E's(i/2-) 


= 1727(06)(06)(11) MeV, 


^=(1/2-) 


= 1825(6)(5)(12) MeV, 



which we present in physical units. Since the lattice spacing is known with less precision 
than the lattice masses presented here, we make the systematic uncertainty arising from the 
lattice spacing explicit (third uncertainty). Given that there are significant ambiguities in 
scale setting, the most precise result will be for dimensionless quantities. 
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With high precision measurements of baryon correlation functions obtained from a single 
type of source for the light-quark and strange-quarks propagators, we have shown that 
the number of methods that can be used to extract the arguments of the contributing 
exponentials increases. This is due to the fact that some methods become stable when 
the uncertainties become small, such as the method of Prony and also the direct fitting of 
multiple exponentials. Histograms constructed from the roots found in the Prony method 
are found to be useful in identifying mass regions where states may exist, but we have not 
yet arrived at a well-defined (rigorous) method with which to use these histograms directly. 
It is likely that a more refined statistical analysis of these correlation functions using the 
most modern statistical tools will further increase the physics that can be extracted. 

The exponential signal-to-noise degradation that plagues baryon correlation functions 
is currently a serious limitation for the calculation of nuclear physics observables, and is 
one significant difference between particle and nuclear physics LQCD calculations. The 
high statistics calculations we have performed have allowed us to systematically explore 
this issue. We find that the issue is more serious than one would naively expect, due 
to (what in hindsight is now obvious) the use of anti-periodic EC's in the time-direction 
on the quark-propagators. The variance of the correlation function is symmetric about 
the mid-point of the time-direction of the configuration. Therefore, the optimal region in 
which to determine the baryon masses (and also their interactions) is in the first half of 
the configuration, far from the midpoint. This significantly reduces the number of useful 
time- slices. Given that most the time required for these calculations is in the measurements, 
and not in the configuration generation, a cure for this problem is to generate ensembles 
of configurations that are longer in the time-direction than those currently being used (as 
opposed to working with different EC's on the quark propagators that arc less theoretically 
"clean" ).^ The multi-exponential fitting and Prony methods enable the ground-state to be 
probed closer to the source where the statistical uncertainties are exponentially smaller (also 
one of the important aspects of the variational method), somewhat reducing the impact of 
the exponentially degrading signal-to-noise near the mid-point of the configuration. Given 
that the signal-to-noisc degradation is exponentially more severe in systems containing two 
or more baryons, all currently available tools will be required to make optimal use of the 
computational resources. For excited states in a given channel, variational methods seem to 
be superior to the standard approach used here. 

An important result of this work has been to quantify the statistical scaling of simple 
observables in the sub-percent regime of uncertainty. Scaling with the number of configu- 
rations was found to adhere to the expected 1 / -^/Ndg behaviour. We have also investigated 
the issue of saturation, asking many measurements can be performed on a single gauge-field 
configuration before it becomes more cost effective to generate another statistically indepen- 
dent gauge- field configuration? To address this we have looked for deviations from 1 / y/N^ 
behavior in the uncertainties in the correlation functions and extracted masses as a function 
of the number of measurements performed on each configuration in the ensemble. The mea- 
surements of the mesons start to saturate after a relatively small number of measurements, 
in this case of order ~ 10, while the baryon correlation functions show no signs of saturation 



^ We note that combining quark propagators with both periodic and anti-periodic temporal BCs, to effec- 
tively double the length of the configurations as seen by the valence quarks, will not resolve all the noise 
issues as much of the problem is produced by states involving sea quarks which are encoded in the gauge 
configurations. 
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up to ~ 200. 

A natural question to ask is if the same extractions could have been performed with 
fewer computational resources by using the variational method with a number of different 
sources for the baryons. We estimate that comparable resources would have been required 
to achieve comparable uncertainties in the states we have examined. However, we have not 
been able to extract excited states of the nucleon with much precision because of closely 
spaced states with the same quantum numbers. This is likely to be something that the 
variational method would better control. Given that the present work was exploratory in 
nature, this is not a concern at present, but it is clear that high-statistics calculations of 
correlation functions arising from multiple interpolating operators will be required in order 
to explore the structure and interactions of nuclei. We are working on implementing this, but 
it will require significant computational resources to perform, even at pion mass ~ 390 MeV 
and for a relatively small lattice volume and relatively coarse lattice spacing. 

It is clear that sub-MeV uncertainties an hadron energies will become routine with the 
anticipated increase in computational resources available to lattice QCD, and that the small 
energy scales that characterize nuclear physics are within reach. However, this program 
will require large ensembles of gauge-field configurations that have large extent in the time- 
direction, and will require a large fraction of the computational resources devoted to mea- 
surements. 
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